home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / fractalobject.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  32.7 KB  |  917 lines
  1. /*
  2.  * fractalobject.c
  3.  *
  4.  * Peter Janssen 
  5.  *
  6.  */
  7. #include "geom.h"
  8. #include "fractalobject.h"
  9. #include "bounds.h"
  10. #include <math.h>
  11. #include <stdio.h>
  12. #include <stdlib.h>
  13.  
  14.  
  15. FractalClosure   *Fractalclosure;
  16. PointPool        *Pointpool;
  17. TrianglePool     *Trianglepool;
  18. EntityPool       *Entitypool;
  19.  
  20.  
  21. static Methods *iFractalObjMethods = NULL;
  22. static char fractalobjName[] = "fractalobject";
  23. static unsigned long raynumber = 1;
  24. static FirstPoint  *SideList;
  25. unsigned long FractalObjTests, FractalObjHits;
  26. static Float Scale;
  27.  
  28.  
  29. /* 
  30.  * Random Generator 
  31.  */
  32.  
  33. #define rand_m  (unsigned long)2147483647
  34. #define rand_q  (unsigned long)127773
  35.  
  36. #define rand_a (unsigned int)16807
  37. #define rand_r (unsigned int)2836
  38.  
  39. /*
  40.  * F(z) = (az)%m
  41.  *      = az-m(az/m)
  42.  *
  43.  * F(z)  = G(z)+mT(z)
  44.  * G(z)  = a(z%q)- r(z/q)
  45.  * T(z)  = (z/q) - (az/m)
  46.  *
  47.  * F(z)  = a(z%q)- rz/q+ m((z/q) - a(z/m))
  48.  *       = a(z%q)- rz/q+ m(z/q) - az
  49.  */
  50.  
  51. static long rand_seed;
  52.  
  53. static Float rand_number()
  54. {
  55. long    lo, hi, test;
  56. static Float inv_rand_m;
  57.     inv_rand_m = 1.0 / rand_m;
  58.  
  59.     hi   = rand_seed/rand_q;
  60.     lo   = rand_seed%rand_q;
  61.  
  62.     test = rand_a*lo - rand_r*hi;
  63.  
  64.     if (test > 0)
  65.         rand_seed = test;
  66.     else
  67.         rand_seed = test + rand_m;
  68.  
  69.     return (rand_seed * inv_rand_m * 2) - 1;
  70. }
  71.  
  72. /*
  73.  * Calculate the boundingbox and the size of a voxel
  74.  * for every subobject in the fractalobject
  75.  */
  76. static void CalculateBoundingBoxesAndVoxelSizes(fractalobj)
  77. FractalObj  *fractalobj;
  78. {
  79.    int            i, j, CurrPointNumber;
  80.    SubObject     *CurrSubobject;
  81.    FractalPoint  *Points;
  82.  
  83.    Points = Fractalclosure->Pointpool->Points;
  84.    for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) {
  85.        for (j = 0; j <= 2; j++) CurrSubobject->Bounds[LOW][j] = MAXVALUE;
  86.        for (j = 0; j <= 2; j++) CurrSubobject->Bounds[HIGH][j] = -MAXVALUE;
  87.        for (i = 0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) {
  88.             for (j = 0; j <= 2; j++) {
  89.                  CurrPointNumber = CurrSubobject->Trianglepool->Triangles[i].Point[j];
  90.                  CurrSubobject->Bounds[LOW][X] = min(CurrSubobject->Bounds[LOW][X], Points[CurrPointNumber].x);
  91.                  CurrSubobject->Bounds[HIGH][X] = max(CurrSubobject->Bounds[HIGH][X], Points[CurrPointNumber].x);
  92.                  CurrSubobject->Bounds[LOW][Y] = min(CurrSubobject->Bounds[LOW][Y], Points[CurrPointNumber].y);
  93.                  CurrSubobject->Bounds[HIGH][Y] = max(CurrSubobject->Bounds[HIGH][Y], Points[CurrPointNumber].y);
  94.                  CurrSubobject->Bounds[LOW][Z] = min(CurrSubobject->Bounds[LOW][Z], Points[CurrPointNumber].z);
  95.                  CurrSubobject->Bounds[HIGH][Z] = max(CurrSubobject->Bounds[HIGH][Z], Points[CurrPointNumber].z);
  96.             }
  97.        }
  98.        CurrSubobject->VoxelSize[X] = (CurrSubobject->Bounds[HIGH][X] - CurrSubobject->Bounds[LOW][X]) /
  99.                                       CurrSubobject->XSize;
  100.        CurrSubobject->VoxelSize[Y] = (CurrSubobject->Bounds[HIGH][Y] - CurrSubobject->Bounds[LOW][Y]) /
  101.                                       CurrSubobject->YSize;
  102.        CurrSubobject->VoxelSize[Z] = (CurrSubobject->Bounds[HIGH][Z] - CurrSubobject->Bounds[LOW][Z]) /
  103.                                       CurrSubobject->ZSize;
  104.    }
  105. }
  106.  
  107.  
  108. /*
  109.  * Calculate new points on the midpoints of the
  110.  * three edges of the triangle
  111.  */
  112. static void CalculatePoints(fractalobj, NewPoints, FTriangle, Delta)
  113. FractalObj       *fractalobj;
  114. int               NewPoints[3];
  115. FractalTriangle  *FTriangle;
  116. DeltaInfo         Delta[3];
  117. {
  118.    int           i, j, fpoint, lpoint, NPoints, Pointsallocated;
  119.    FirstPoint   *Currfp, *chasingfp, *Newfp;
  120.    LastPoint    *Currlp, *chasinglp, *Newlp;
  121.    FractalPoint *Points, *TPoints;
  122.    PointPool    *Pointpool;
  123.    Float         displacement, size1, size2;
  124.    Vector        n1, n2;
  125.  
  126.    Pointpool = Fractalclosure->Pointpool;
  127.    NPoints = Pointpool->NumberOfPoints;
  128.    Pointsallocated = Pointpool->PointsAllocated;
  129.    Points = Pointpool->Points;
  130.    for(i=0; i<=2; i++) {
  131.        fpoint = min(FTriangle->Point[(i + 1)%3], FTriangle->Point[i]);
  132.        lpoint = max(FTriangle->Point[(i + 1)%3], FTriangle->Point[i]);
  133.        for(Currfp = SideList; fpoint > Currfp->Point; chasingfp = Currfp,
  134.               Currfp = Currfp->Next);
  135.        if (Currfp->Point == fpoint) {
  136.            /* 
  137.             * sidelist has an entry with firstpoint 
  138.             */
  139.            chasinglp = NULL;
  140.            for(Currlp = Currfp->Info; Currlp && (Currlp->Point < lpoint);
  141.                   chasinglp = Currlp, Currlp = Currlp->Next);
  142.            if (Currlp && (Currlp->Point == lpoint)) {
  143.                     /* 
  144.                      * sidelist has an entry for pointpair fpoint, lpoint 
  145.                      */
  146.                     NewPoints[i] = Currlp->NewPoint;
  147.                     continue;
  148.            }
  149.        }
  150.  
  151.        /*
  152.         * Calculate New Point 
  153.         */
  154.        displacement = Delta[i].distance * pow((fractalobj->Excentricity * Delta[i].distance * Scale), fractalobj->Dimension);
  155.        /* n1 */
  156.        if (Delta[i].z != 0) {
  157.            n1.x = -Delta[i].z;
  158.            n1.y = -Delta[i].z;
  159.            n1.z =  Delta[i].x + Delta[i].y;
  160.        } else if (Delta[i].y != 0) {
  161.            n1.x = -Delta[i].y;
  162.            n1.y =  Delta[i].x + Delta[i].z;
  163.            n1.z = -Delta[i].y;
  164.        } else {
  165.            n1.x =  Delta[i].y + Delta[i].z;
  166.            n1.y = -Delta[i].x;
  167.            n1.z = -Delta[i].x;
  168.        }
  169.        VecNormalize(&n1);
  170.        /* n2 */
  171.        VecCross(&Delta[i], &n1, &n2);
  172.        VecNormalize(&n2);
  173.        size1 = rand_number();
  174.        size2 = rand_number();
  175.        CheckPoints(Points, NPoints, Pointsallocated);
  176.        Points[NPoints].x = (Points[fpoint].x + Points[lpoint].x) / 2 +
  177.                             displacement * size1 * n1.x + displacement * size2 * n2.x;
  178.        Points[NPoints].y = (Points[fpoint].y + Points[lpoint].y) / 2 +
  179.                             displacement * size1 * n1.y + displacement * size2 * n2.y;
  180.        Points[NPoints].z = (Points[fpoint].z + Points[lpoint].z) / 2 +
  181.                             displacement * size1 * n1.z + displacement * size2 * n2.z;
  182.        NewPoints[i] = NPoints;
  183.  
  184.        /* 
  185.         * Currfp points to position equal to or just bigger than fpoint 
  186.         */
  187.        Newlp = (LastPoint *)share_malloc(sizeof(LastPoint));
  188.        Newlp->Point = lpoint;
  189.        if (Currfp->Point == fpoint) {
  190.                   /* 
  191.                    * Currlp points tot position equal to or just bigger than
  192.                    * lpoint or is equal to NULL 
  193.                    * Insert between chasinglp and Currlp 
  194.                    */
  195.                   if (chasinglp == NULL)
  196.                        /* 
  197.                         * the first in line 
  198.                         */
  199.                        Currfp->Info = Newlp;
  200.                   else
  201.                        chasinglp->Next = Newlp;
  202.                   Newlp->Next = Currlp;
  203.         } else {
  204.             Newfp = (FirstPoint *)share_malloc(sizeof(FirstPoint));
  205.             chasingfp->Next = Newfp;
  206.             Newfp->Next     = Currfp;
  207.             Newfp->Info     = Newlp;
  208.             Newfp->Point    = fpoint;
  209.             Newlp->Next     = NULL;
  210.         }
  211.         Newlp->NewPoint = NPoints++;
  212.    }
  213.    Pointpool->NumberOfPoints  = NPoints;
  214.    Pointpool->PointsAllocated = Pointsallocated;
  215.    Pointpool->Points          = Points;
  216. }
  217.  
  218.  
  219. /* 
  220.  * Take triangulated object and distort it
  221.  */
  222. static void BuildFractalObject(fractalobj)
  223. FractalObj    *fractalobj;
  224. {
  225.     Float             Maxsize, Delta_x, Delta_y, Delta_z, distance;
  226.     int               Changed, Poolchanged, ONTriangles, NTriangles, Trianglesallocated;
  227.     FractalPoint      *FPoints;
  228.     SubObject         *CurrSubobject;
  229.     TrianglePool      *Trianglepool;
  230.     FractalTriangle   *FTriangles;
  231.     DeltaInfo          Delta[3];
  232.     int                i, j, k, NewPoints[3], t1, t2;
  233.     FirstPoint        *Currfp, *chasingfp;
  234.     LastPoint         *Currlp, *chasinglp;
  235.  
  236.     Maxsize = fractalobj->MaximumSize;
  237.     Changed = TRUE;
  238.     Scale = 0;
  239.     FPoints = Pointpool->Points;
  240.     for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) {
  241.              FTriangles = CurrSubobject->Trianglepool->Triangles;
  242.              for (i=0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) {
  243.                   for(j = 0; j <= 2; j++) {
  244.                       t1 = FTriangles[i].Point[j];
  245.                       t2 = FTriangles[i].Point[(j+1)%3];
  246.                       Delta_x = FPoints[t1].x - FPoints[t2].x;
  247.                       Delta_y = FPoints[t1].y - FPoints[t2].y;
  248.                       Delta_z = FPoints[t1].z - FPoints[t2].z;
  249.                       distance = sqrt(Delta_x * Delta_x +
  250.                                       Delta_y * Delta_y +
  251.                                       Delta_z * Delta_z);
  252.                       Scale = max(Scale, distance);
  253.                   } 
  254.              }
  255.     }
  256.     Scale = 1.0 / Scale;
  257.  
  258.     while (Changed) {
  259.         Changed = FALSE;
  260.         for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) {
  261.              Trianglepool = CurrSubobject->Trianglepool;
  262.              if (Trianglepool->PoolChanged) {
  263.                  Poolchanged = FALSE;
  264.                  ONTriangles = NTriangles = Trianglepool->NumberOfTriangles;
  265.                  Trianglesallocated = Trianglepool->TrianglesAllocated;
  266.                  FTriangles = Trianglepool->Triangles;
  267.                  for (i=0; i < ONTriangles; i++) {
  268.                       FPoints = Pointpool->Points;
  269.                       for(j = 0; j <= 2; j++) {
  270.                           t1 = FTriangles[i].Point[j];
  271.                           t2 = FTriangles[i].Point[(j+1)%3];
  272.                           Delta[j].x = FPoints[t1].x - FPoints[t2].x;
  273.                           Delta[j].y = FPoints[t1].y - FPoints[t2].y;
  274.                           Delta[j].z = FPoints[t1].z - FPoints[t2].z;
  275.                           Delta[j].distance = sqrt(Delta[j].x * Delta[j].x +
  276.                                                    Delta[j].y * Delta[j].y +
  277.                                                    Delta[j].z * Delta[j].z);
  278.                       }
  279.                       if ((Delta[0].distance > Maxsize) ||
  280.                           (Delta[1].distance > Maxsize) ||
  281.                           (Delta[2].distance > Maxsize)) {
  282.                              Poolchanged = Changed = TRUE;
  283.                              if (NTriangles + 3 >= Trianglesallocated) {
  284.                                 Trianglesallocated += TrianglesToAllocate;
  285.                                 FTriangles = (FractalTriangle *)realloc((char *)FTriangles, Trianglesallocated * sizeof(FractalTriangle));
  286.                              }
  287.                              NTriangles += 3;
  288.                              CalculatePoints(fractalobj, NewPoints, &FTriangles[i], Delta);
  289.                              FTriangles[NTriangles - 1].Point[0] = NewPoints[0];
  290.                              FTriangles[NTriangles - 1].Point[1] = NewPoints[1];
  291.                              FTriangles[NTriangles - 1].Point[2] = NewPoints[2];
  292.                              FTriangles[NTriangles - 3].Point[0] = NewPoints[0];
  293.                              FTriangles[NTriangles - 3].Point[1] = FTriangles[i].Point[1];
  294.                              FTriangles[NTriangles - 3].Point[2] = NewPoints[1];
  295.                              FTriangles[NTriangles - 2].Point[0] = NewPoints[2];
  296.                              FTriangles[NTriangles - 2].Point[1] = NewPoints[1];
  297.                              FTriangles[NTriangles - 2].Point[2] = FTriangles[i].Point[2];
  298.                              FTriangles[i].Point[1] = NewPoints[0];
  299.                              FTriangles[i].Point[2] = NewPoints[2];
  300.                       }
  301.                  }
  302.                  Trianglepool->PoolChanged        = Poolchanged;
  303.                  Trianglepool->NumberOfTriangles  = NTriangles;
  304.                  Trianglepool->TrianglesAllocated = Trianglesallocated;
  305.                  Trianglepool->Triangles          = FTriangles;
  306.             }
  307.         }
  308.         /* 
  309.          * release sidelist and setup new one 
  310.          */
  311.         for(Currfp = SideList->Next; Currfp->Point != MAXVALUE; chasingfp = Currfp, Currfp = Currfp->Next, free(chasingfp))
  312.             for(Currlp = Currfp->Info; Currlp; chasinglp = Currlp,
  313.                      Currlp = Currlp->Next, free(chasinglp));
  314.         SideList->Next = Currfp;
  315.     }
  316.     CalculateBoundingBoxesAndVoxelSizes(fractalobj);
  317. }
  318.  
  319.  
  320. /*
  321.  * Convert worldcoordinates to gridcoordinates 
  322.  */
  323. static void Pos2Grid(subobj, pos, index)
  324. SubObject   *subobj;
  325. Float        pos[3];
  326. int          index[3];
  327. {
  328.  
  329.         index[X] = (int)(Fx2voxel(subobj, pos[0]));
  330.         index[Y] = (int)(Fy2voxel(subobj, pos[1]));
  331.         index[Z] = (int)(Fz2voxel(subobj, pos[2]));
  332.  
  333.         if (index[X] == subobj->XSize)
  334.                 index[X]--;
  335.         if (index[Y] == subobj->YSize)
  336.                 index[Y]--;
  337.         if (index[Z] == subobj->ZSize)
  338.                 index[Z]--;
  339. }
  340.  
  341.  
  342.  
  343. /*
  344.  * Put triangles of fractalobject in voxels
  345.  * for efficient raytracing
  346.  */
  347. static void ConvertFractalObject(fractalobj)
  348. FractalObj   *fractalobj;
  349. {
  350.    int               i, x, y, z, low[3], high[3];
  351.    Float             bounds[2][3];
  352.    SubObject        *CurrSubobject;
  353.    FractalTriangle  *Triangles, *CurrTriangle;
  354.    TriangleInfo     *Triangleinfo;
  355.    TriangleClosure  *Triangleclosure;
  356.    FractalPoint     *Points;
  357.    int              p0, p1, p2;
  358.  
  359.    Points = Fractalclosure->Pointpool->Points;
  360.    for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) {
  361.        Triangles = CurrSubobject->Trianglepool->Triangles;
  362.        CurrSubobject->Cells = (TriangleClosure ****)share_malloc(CurrSubobject->XSize *
  363.                     sizeof(TriangleClosure ***));
  364.        for (x = 0; x < CurrSubobject->XSize; x++) {
  365.              CurrSubobject->Cells[x] = (TriangleClosure ***)share_malloc(CurrSubobject->YSize *
  366.                 sizeof(TriangleClosure **));
  367.              for (y = 0; y < CurrSubobject->YSize; y++) 
  368.                   CurrSubobject->Cells[x][y] = (TriangleClosure **)share_calloc(
  369.                     CurrSubobject->ZSize , sizeof(TriangleClosure *));
  370.        }
  371.        for (i = 0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) {
  372.             CurrTriangle = &Triangles[i];
  373.             p0 = CurrTriangle->Point[0];
  374.             p1 = CurrTriangle->Point[1];
  375.             p2 = CurrTriangle->Point[2];
  376.             bounds[LOW][X] = min3(Points[p0].x, Points[p1].x, Points[p2].x);
  377.             bounds[HIGH][X] = max3(Points[p0].x, Points[p1].x, Points[p2].x);
  378.             bounds[LOW][Y] = min3(Points[p0].y, Points[p1].y, Points[p2].y);
  379.             bounds[HIGH][Y] = max3(Points[p0].y, Points[p1].y, Points[p2].y);
  380.             bounds[LOW][Z] = min3(Points[p0].z, Points[p1].z, Points[p2].z);
  381.             bounds[HIGH][Z] = max3(Points[p0].z, Points[p1].z, Points[p2].z);
  382.  
  383.             Pos2Grid(CurrSubobject, bounds[LOW], low);
  384.             Pos2Grid(CurrSubobject, bounds[HIGH], high);
  385.  
  386.             Triangleinfo = (TriangleInfo *)share_malloc(sizeof(TriangleInfo));
  387.             Triangleinfo->Point[0] = p0;
  388.             Triangleinfo->Point[1] = p1;
  389.             Triangleinfo->Point[2] = p2;
  390.             Triangleinfo->Counter = 0;
  391.             for (x = low[X]; x <= high[X]; x++)
  392.                  for (y = low[Y]; y <= high[Y]; y++)
  393.                       for (z = low[Z]; z <= high[Z]; z++) {
  394.                            Triangleclosure = (TriangleClosure *)share_malloc(sizeof(TriangleClosure));
  395.                            Triangleclosure->Triangle = Triangleinfo;
  396.                            Triangleclosure->Next = CurrSubobject->Cells[x][y][z];
  397.                            CurrSubobject->Cells[x][y][z] = Triangleclosure;
  398.                       }
  399.        }
  400.        /* 
  401.         * release trianglepool
  402.         */
  403.        free(CurrSubobject->Trianglepool->Triangles);
  404.        free(CurrSubobject->Trianglepool);
  405.    }
  406. }
  407.  
  408.  
  409. /*
  410.  * Create FractalObject and return reference to it
  411.  */
  412. FractalObj *FractalObjCreate(Number, Dimension, MaxSize, Excentricity, fc)
  413. Float            MaxSize, Dimension, Excentricity;
  414. int              Number;
  415. FractalClosure  *fc;
  416. {
  417.         FractalObj       *fractalobj;
  418.         int              *TList, TriangleNumber;
  419.         FractalTriangle  *InitialTriangles, *LocalTriangles, *TTriangles;
  420.         FractalEntity    *CurrEntity;
  421.         int               NPoints, NTriangles, NEntities, i, j, LocalAllocated, LocalNTriangles;
  422.         SubObject        *Subobjects, *TSubobject;
  423.         TriangleList     *CurrTriangle, *chasingTriangle;
  424.         int               CurrPoint;
  425.  
  426.  
  427.         fractalobj = (FractalObj *)share_malloc(sizeof(FractalObj));
  428.         fractalobj->MaximumSize = MaxSize;
  429.         fractalobj->Dimension = Dimension;
  430.         fractalobj->Excentricity = Excentricity;
  431.         fractalobj->Closure = fc;
  432.         rand_seed = Number;
  433.  
  434.         Fractalclosure = fc;
  435.  
  436.         Pointpool    = Fractalclosure->Pointpool;
  437.         Trianglepool = Fractalclosure->Trianglepool;
  438.         Entitypool   = Fractalclosure->Entitypool;
  439.         NPoints      = Pointpool->NumberOfPoints;
  440.         NTriangles   = Trianglepool->NumberOfTriangles;
  441.         NEntities    = Entitypool->NumberOfEntities;
  442.  
  443.         if (NPoints <= 0)
  444.              RLerror(RL_ABORT, "Illegal number of Points : %d", NPoints);
  445.  
  446.         if (NTriangles <= 0)
  447.              RLerror(RL_ABORT, "Illegal number of Triangles : %d", NTriangles);
  448.  
  449.         if (NEntities <= 0)
  450.              RLerror(RL_ABORT, "Illegal number of Entities : %d", NEntities);
  451.  
  452.         TList = (int *)Malloc(sizeof(int) * NTriangles);
  453.  
  454.         Subobjects = NULL;
  455.         InitialTriangles = Trianglepool->Triangles;
  456.         for (i=0; i < NEntities; i++) {
  457.             CurrEntity = &Entitypool->Entities[i];
  458.             TSubobject = (SubObject *)share_malloc(sizeof(SubObject));
  459.             TSubobject->Trianglepool = (TrianglePool *)share_malloc(sizeof(TrianglePool));
  460.             LocalTriangles = (FractalTriangle *)share_malloc(TrianglesToAllocate * sizeof(FractalTriangle));
  461.             LocalAllocated = TrianglesToAllocate;
  462.             TSubobject->XSize = CurrEntity->XSize;
  463.             TSubobject->YSize = CurrEntity->YSize;
  464.             TSubobject->ZSize = CurrEntity->ZSize;
  465.             LocalNTriangles = 0;
  466.             for (CurrTriangle = CurrEntity->Triangles; CurrTriangle;
  467.                    chasingTriangle = CurrTriangle, CurrTriangle = CurrTriangle->Next,
  468.                    free(chasingTriangle)) {
  469.                  TriangleNumber = CurrTriangle->TriangleNumber;
  470.  
  471.                  if (TriangleNumber < 0 || TriangleNumber >= NTriangles)
  472.                         RLerror(RL_ABORT, "Illegal Triangle number : %d\n", TriangleNumber+1);
  473.                  CheckTriangles(LocalTriangles, LocalNTriangles, LocalAllocated);
  474.  
  475.                  CurrPoint = InitialTriangles[TriangleNumber].Point[0];
  476.                  if ((CurrPoint < 0) || (CurrPoint >= NPoints))
  477.                         RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1);
  478.                  LocalTriangles[LocalNTriangles].Point[0] = CurrPoint;
  479.  
  480.                  CurrPoint = InitialTriangles[TriangleNumber].Point[1];
  481.                  if ((CurrPoint < 0) || (CurrPoint >= NPoints))
  482.                         RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1);
  483.                  LocalTriangles[LocalNTriangles].Point[1] = CurrPoint;
  484.  
  485.                  CurrPoint = InitialTriangles[TriangleNumber].Point[2];
  486.                  if ((CurrPoint < 0) || (CurrPoint >= NPoints))
  487.                         RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1);
  488.                  LocalTriangles[LocalNTriangles].Point[2] = CurrPoint;
  489.  
  490.                  LocalNTriangles++;
  491.                  TList[TriangleNumber] = TRUE;
  492.             }
  493.             TSubobject->Trianglepool->NumberOfTriangles = LocalNTriangles;
  494.             TSubobject->Trianglepool->TrianglesAllocated = LocalAllocated;
  495.             TSubobject->Trianglepool->Triangles = LocalTriangles;
  496.             TSubobject->Trianglepool->PoolChanged = TRUE;
  497.             TSubobject->Next = Subobjects;
  498.             Subobjects = TSubobject;
  499.         }
  500.         fractalobj->SubObjects = Subobjects;
  501.         /* 
  502.          * Check if every triangle is placed in a subobject 
  503.          */
  504.         for (i=0; i < NTriangles; i++)
  505.              if (!TList[i]) RLerror(RL_ADVISE, "Triangle %d not placed in entity.\n", i+1);
  506.         free(TList);
  507.         /* 
  508.          * release triangles uit trianglepool 
  509.          */
  510.         free(Trianglepool->Triangles);
  511.         free(Trianglepool);
  512.         /* 
  513.          * release entitypool 
  514.          */
  515.         free(Entitypool->Entities);
  516.         free(Entitypool);
  517.         /* 
  518.          * Build empty sidelist 
  519.          */
  520.         SideList = (FirstPoint *)share_malloc(sizeof(FirstPoint));
  521.         SideList->Point = -1;
  522.         SideList->Next = (FirstPoint *)share_malloc(sizeof(FirstPoint));
  523.         SideList->Next->Point = MAXVALUE;
  524.  
  525.         /* 
  526.          * Distort object 
  527.          */
  528.         BuildFractalObject(fractalobj);
  529.         /* 
  530.          * free sidelist
  531.          */
  532.         free(SideList->Next);
  533.         free(SideList);
  534.         /* 
  535.          * Convert to voxels 
  536.          */
  537.         ConvertFractalObject(fractalobj);
  538.  
  539.         return fractalobj;
  540. }
  541.  
  542.  
  543.  
  544. /*
  545.  * Intersect ray with fractalobject.
  546.  */
  547. int FractalObjIntersect(fractalobj, ray, mindist, maxdist)
  548. FractalObj *fractalobj;
  549. Ray *ray;
  550. Float mindist, *maxdist;
  551. {
  552.         int            hit;
  553.         Float          offset;
  554.         Vector         curpos;
  555.         unsigned long  counter;
  556.         SubObject     *CurrSubobject;
  557.  
  558.         FractalObjTests++;
  559.         hit = FALSE;
  560.  
  561.         /*
  562.          * If outside of the bounding box, check to see if we
  563.          * hit it.
  564.          */
  565.         VecAddScaled(ray->pos, mindist, ray->dir, &curpos);
  566.         offset = *maxdist;
  567.         if (OutOfBounds(&curpos, fractalobj->Bounds) &&
  568.             !BoundsIntersect(ray, fractalobj->Bounds, mindist, &offset))
  569.                         /*
  570.                          * Ray never hit grid space.
  571.                          */
  572.                         return hit;
  573.  
  574.         counter = raynumber++;
  575.  
  576.         for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next)
  577.              if (SubObjectIntersect(fractalobj, CurrSubobject, ray, mindist, maxdist, counter))
  578.                     hit = TRUE;
  579.        return hit;
  580. }
  581.  
  582.  
  583. /*
  584.  * Test for intersection of ray with every of the subobjects
  585.  *
  586.  * (based on the code of grid.c)
  587.  */
  588. int SubObjectIntersect(fractalobj, subobject, ray, mindist, maxdist, counter)
  589. FractalObj      *fractalobj;
  590. SubObject       *subobject;
  591. Ray             *ray;
  592. Float            mindist;
  593. Float           *maxdist;
  594. unsigned long    counter;
  595. {
  596.         int hit;
  597.         Float offset, tMaxX, tMaxY, tMaxZ;
  598.         Float tDeltaX, tDeltaY, tDeltaZ;
  599.         int stepX, stepY, stepZ, outX, outY, outZ, x, y, z;
  600.         Vector curpos;
  601.         TriangleClosure  *list;
  602.  
  603.         hit = FALSE;
  604.  
  605.         /*
  606.          * If outside of the bounding box, check to see if we
  607.          * hit it.
  608.          */
  609.         VecAddScaled(ray->pos, mindist, ray->dir, &curpos);
  610.         if (OutOfBounds(&curpos, subobject->Bounds)) {
  611.                 offset = *maxdist;
  612.                 if (!BoundsIntersect(ray, subobject->Bounds, mindist, &offset))
  613.                         /*
  614.                          * Ray never hit grid space.
  615.                          */
  616.                         return hit;
  617.                 VecAddScaled(ray->pos, offset, ray->dir, &curpos);
  618.         } else
  619.                 offset = mindist;
  620.  
  621.  
  622.         /*
  623.          * tMaxX is the absolute distance from the ray origin we must move
  624.          *              to get to the next voxel in the X
  625.          *              direction.  It is incrementally updated
  626.          *              by DDA as we move from voxel-to-voxel.
  627.          * tDeltaX is the relative amount along the ray we move to
  628.          *              get to the next voxel in the X direction. Thus,
  629.          *              when we decide to move in the X direction,
  630.          *              we increment tMaxX by tDeltaX.
  631.          */
  632.         x = Fx2voxel(subobject, curpos.x);
  633.         if (x == subobject->XSize)
  634.                 x--;
  635.         if (fabs(ray->dir.x) < EPSILON) {
  636.                 tMaxX = FAR_AWAY;
  637.                 tDeltaX = 0.;
  638.         } else if (ray->dir.x < 0.) {
  639.                 tMaxX = offset + (Fvoxel2x(subobject, x) - curpos.x) / ray->dir.x;
  640.                 tDeltaX = subobject->VoxelSize[X] / - ray->dir.x;
  641.                 stepX = outX = -1;
  642.         } else {
  643.                 tMaxX = offset + (Fvoxel2x(subobject, x+1) - curpos.x) / ray->dir.x;
  644.                 tDeltaX = subobject->VoxelSize[X] / ray->dir.x;
  645.                 stepX = 1;
  646.                 outX = subobject->XSize;
  647.         }
  648.  
  649.         y = Fy2voxel(subobject, curpos.y);
  650.         if (y == subobject->YSize)
  651.                 y--;
  652.  
  653.         if (fabs(ray->dir.y) < EPSILON) {
  654.                 tMaxY = FAR_AWAY;
  655.                 tDeltaY = 0.;
  656.         } else if (ray->dir.y < 0.) {
  657.                 tMaxY = offset + (Fvoxel2y(subobject, y) - curpos.y) / ray->dir.y;
  658.                 tDeltaY = subobject->VoxelSize[Y] / - ray->dir.y;
  659.                 stepY = outY = -1;
  660.         } else {
  661.                 tMaxY = offset + (Fvoxel2y(subobject, y+1) - curpos.y) / ray->dir.y;
  662.                 tDeltaY = subobject->VoxelSize[Y] / ray->dir.y;
  663.                 stepY = 1;
  664.                 outY = subobject->YSize;
  665.         }
  666.  
  667.         z = Fz2voxel(subobject, curpos.z);
  668.         if (z == subobject->ZSize)
  669.                 z--;
  670.         if (fabs(ray->dir.z) < EPSILON) {
  671.                 tMaxZ = FAR_AWAY;
  672.                 tDeltaZ = 0.;
  673.         } else if (ray->dir.z < 0.) {
  674.                 tMaxZ = offset + (Fvoxel2z(subobject, z) - curpos.z) / ray->dir.z;
  675.                 tDeltaZ = subobject->VoxelSize[Z] / - ray->dir.z;
  676.                 stepZ = outZ = -1;
  677.         } else {
  678.                 tMaxZ = offset + (Fvoxel2z(subobject, z+1) - curpos.z) / ray->dir.z;
  679.                 tDeltaZ = subobject->VoxelSize[Z] / ray->dir.z;
  680.                 stepZ = 1;
  681.                 outZ = subobject->ZSize;
  682.         }
  683.  
  684.         /* 
  685.          * Do DDA-algoritm to traverse the grid 
  686.          */
  687.         while (TRUE) {
  688.                 list = subobject->Cells[x][y][z];
  689.                 if (tMaxX < tMaxY && tMaxX < tMaxZ) {
  690.                         if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist))
  691.                                         hit = TRUE;
  692.                         x += stepX;
  693.                         if (*maxdist < tMaxX || x == outX)
  694.                                 break;
  695.                         tMaxX += tDeltaX;
  696.                 } else if (tMaxZ < tMaxY) {
  697.                         if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist))
  698.                                         hit = TRUE;
  699.                         z += stepZ;
  700.                         if (*maxdist < tMaxZ || z == outZ)
  701.                                 break;
  702.                         tMaxZ += tDeltaZ;
  703.                 } else {
  704.                         if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist))
  705.                                         hit = TRUE;
  706.                         y += stepY;
  707.                         if (*maxdist < tMaxY || y == outY)
  708.                                 break;
  709.                         tMaxY += tDeltaY;
  710.                 }
  711.         }
  712.         return hit;
  713. }
  714.  
  715. /*
  716.  * Do intersection test with every triangle in the voxel 
  717.  */
  718. static int ListIntersect(fractalobj, list, ray, counter, mindist, maxdist)
  719. FractalObj      *fractalobj;
  720. TriangleClosure *list;
  721. Ray             *ray;
  722. unsigned long    counter;
  723. Float            mindist, *maxdist;
  724. {
  725.         int hit;
  726.  
  727.         hit = FALSE;
  728.  
  729.         for(;list; list = list->Next) {
  730.                 /*
  731.                  * If object's counter is greater than or equal to the
  732.                  * number associated with the current grid,
  733.                  * don't bother checking again.
  734.                  */
  735.                 if (list->Triangle->Counter < counter ) {
  736.                         list->Triangle->Counter = counter;
  737.                         if (TriangleIntersect(fractalobj, list->Triangle->Point, ray, mindist, maxdist))
  738.                                 hit = TRUE;
  739.                 }
  740.         }
  741.         return hit;
  742. }
  743.  
  744.  
  745. /*
  746.  * Do intersection test with triangle 
  747.  */
  748. static int TriangleIntersect(fractalobj, point, ray, mindist, maxdist)
  749. FractalObj *fractalobj;
  750. int         point[3];
  751. Ray        *ray;
  752. Float       mindist, *maxdist;
  753. {
  754.     Vector  edge1, edge2, edge3;
  755.     Float   k, s;
  756.     Vector  normal, absnorm;
  757.     Vector  pos, dir;
  758.     Float   qi1, qi2, b0, b1, b2;
  759.     FractalPoint  *Points;
  760.     FractalPoint  p0, p1, p2;
  761.  
  762.  
  763.     Points = fractalobj->Closure->Pointpool->Points;
  764.     p0 = Points[point[0]];
  765.     p1 = Points[point[1]];
  766.     p2 = Points[point[2]];
  767.     VecSub(p1, p0, &edge1);
  768.     VecSub(p2, p1, &edge2);
  769.     VecSub(p0, p2, &edge3);
  770.  
  771.     /*
  772.      *  Find plane normal. 
  773.      */
  774.     VecCross(&edge1, &edge2, &normal);
  775.  
  776.     pos = ray->pos;
  777.     dir = ray->dir;
  778.     /*
  779.      * Plane intersection.
  780.      */
  781.     k = dotp(&normal, &dir);
  782.     if (fabs(k) < EPSILON)
  783.         return FALSE;
  784.     s = (dotp(&normal, &p0) - dotp(&normal, &pos)) / k;
  785.     if (s < mindist || s > *maxdist)
  786.                 return FALSE;
  787.  
  788.     /*
  789.      * Find "dominant" part of normal vector.
  790.      */
  791.     absnorm.x = fabs(normal.x);
  792.     absnorm.y = fabs(normal.y);
  793.     absnorm.z = fabs(normal.z);
  794.  
  795.     if (absnorm.x > absnorm.y && absnorm.x > absnorm.z) {
  796.         qi1 = pos.y + s * dir.y;
  797.         qi2 = pos.z + s * dir.z;
  798.         b0 = (edge2.y * (qi2 - p1.z) -
  799.               edge2.z * (qi1 - p1.y)) / normal.x;
  800.         if (b0 < 0. || b0 > 1.) return FALSE;
  801.         b1 = (edge3.y * (qi2 - p2.z) -
  802.               edge3.z * (qi1 - p2.y)) / normal.x;
  803.         if (b1 < 0. || b1 > 1.) return FALSE;
  804.         b2 = (edge1.y * (qi2 - p0.z) -
  805.               edge1.z * (qi1 - p0.y)) / normal.x;
  806.         if (b2 < 0. || b2 > 1.) return FALSE;
  807.     } else if (absnorm.y > absnorm.z) {
  808.         qi1 = pos.x + s * dir.x;
  809.         qi2 = pos.z + s * dir.z;
  810.         b0 = (edge2.z * (qi1 - p1.x) -
  811.               edge2.x * (qi2 - p1.z)) / normal.y;
  812.         if (b0 < 0. || b0 > 1.) return FALSE;
  813.         b1 = (edge3.z * (qi1 - p2.x) -
  814.               edge3.x * (qi2 - p2.z)) / normal.y;
  815.         if (b1 < 0. || b1 > 1.) return FALSE;
  816.         b2 = (edge1.z * (qi1 - p0.x) -
  817.               edge1.x * (qi2 - p0.z)) / normal.y;
  818.         if (b2 < 0. || b2 > 1.) return FALSE;
  819.     } else {
  820.         qi1 = pos.x + s * dir.x;
  821.         qi2 = pos.y + s * dir.y;
  822.         b0 = (edge2.x * (qi2 - p1.y) -
  823.               edge2.y * (qi1 - p1.x)) / normal.z;
  824.         if (b0 < 0. || b0 > 1.) return FALSE;
  825.         b1 = (edge3.x * (qi2 - p2.y) -
  826.               edge3.y * (qi1 - p2.x)) / normal.z;
  827.         if (b1 < 0. || b1 > 1.) return FALSE;
  828.         b2 = (edge1.x * (qi2 - p0.y) -
  829.               edge1.y * (qi1 - p0.x)) / normal.z;
  830.         if (b2 < 0. || b2 > 1.) return FALSE;
  831.     }
  832.     FractalObjHits++;
  833.     *maxdist = s;
  834.     fractalobj->Normal = normal;
  835.     VecAddScaled(ray->pos, s, ray->dir, &fractalobj->IntersectPos);
  836.  
  837.     return TRUE;
  838. }
  839.  
  840. #define VecEqual(a,b) ((equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z)))
  841.  
  842. /*
  843.  * Calculate Normal on triangle at pos
  844.  */
  845. int FractalObjNormal(fractalobj, pos, nrm, gnrm)
  846. FractalObj *fractalobj;
  847. Vector     *pos, *nrm, *gnrm;
  848. {
  849.     if (!VecEqual(*pos, fractalobj->IntersectPos))
  850.         RLerror(RL_ABORT, "Can't compute fractalobj normal in point other than the last intersection.\n");
  851.         *gnrm = fractalobj->Normal;
  852.  
  853.         VecNormalize(gnrm);
  854.         *nrm = *gnrm;
  855.         return FALSE;
  856. }
  857.  
  858.  
  859. char *FractalObjName()
  860. {
  861.         return fractalobjName;
  862. }
  863.  
  864. void FractalObjStats(tests, hits)
  865. unsigned long *tests, *hits;
  866. {
  867.         *tests = FractalObjTests;
  868.         *hits = FractalObjHits;
  869. }
  870.  
  871. void FractalObjBounds(fractalobj, bounds)
  872. FractalObj  *fractalobj;
  873. Float        bounds[2][3];
  874. {
  875.     int         i, j;
  876.     SubObject  *CurrSubobject;
  877.  
  878.     for (i = 0; i <= 2; i++) {
  879.          bounds[LOW][i] = MAXVALUE;
  880.          bounds[HIGH][i] = -MAXVALUE;
  881.     }
  882.  
  883.     for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next)
  884.         for (i = 0; i <= 2; i++) {
  885.              bounds[LOW][i] = min(bounds[LOW][i], CurrSubobject->Bounds[LOW][i]);
  886.              bounds[HIGH][i] = max(bounds[HIGH][i], CurrSubobject->Bounds[HIGH][i]);
  887.         }
  888.     for (i = 0; i <= 2; i++)
  889.          for (j = 0; j <= 1; j++)
  890.               fractalobj->Bounds[j][i] = bounds[j][i];
  891. }
  892.  
  893. void FractalObjMethodRegister(meth)
  894. UserMethodType  meth;
  895. {
  896.         if (iFractalObjMethods)
  897.             iFractalObjMethods->user = meth;
  898. }
  899.  
  900.  
  901. Methods *FractalObjMethods()
  902. {
  903.         if (iFractalObjMethods == (Methods *)NULL) {
  904.                 iFractalObjMethods = MethodsCreate();
  905.                 iFractalObjMethods->create = (GeomCreateFunc *)FractalObjCreate;
  906.                 iFractalObjMethods->methods = FractalObjMethods;
  907.                 iFractalObjMethods->name = FractalObjName;
  908.                 iFractalObjMethods->intersect = FractalObjIntersect;
  909.                 iFractalObjMethods->normal = FractalObjNormal;
  910.                 iFractalObjMethods->bounds = FractalObjBounds;
  911.                 iFractalObjMethods->stats = FractalObjStats;
  912.                 iFractalObjMethods->checkbounds = FALSE;
  913.                 iFractalObjMethods->closed = FALSE;
  914.         }
  915.         return iFractalObjMethods;
  916. }
  917.